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ABSTRACT 

We discuss a possible form for a theory akin to local density functional the- 
ory, but able to produce van der Waals energies in a natural fashion. The 
usual Local Density Approximation (LDA) for the exchange and correlation 
energy E^c of an inhomogeneous electronic system can be derived by making 
a quasilocal approximation for the interacting density-density response func- 
tion xC?', then using the fluctuation-dissipation theorem and a Feynman 
coupling-constant integration to generate Exc- The first new idea proposed here 
is to use the same approach except that one makes a quasilocal approximation 
for the bare response rather than for x- The interacting response is then 
obtained by solving a nonlocal screening integral equation in real space. If the 
nonlocal screening is done at the time-dependent Hartree level, then the re- 
sulting energy is an approximation to the full inhomogeneous RPA energy: we 
show here that the inhomogeneous RPA correlation energy contains a van der 
Waals term for the case of widely-separated neutral subsystems. The second 
new idea is to use a particularly simple way of introducing LDA-like local field 
corrrections into the screening equations, giving a theory which should remain 
reasonable for all separations of a pair of subsystems, encompassing both the 
van der Waals limit much as in RPA and the bonding limit much as in LDA 
theory. The resulting functional is an explicit recipe which takes a trial electron 
density distribution n{r) as its input and yields E^c as its output. The reason 
that it includes the van der Waals energy is that, like the RPA formalism but 
unlike LDA, it involves spatially nonlocal screening for inhomogeneous systems. 
It differs from earlier theories which effectively use hydrodnamics to generate 
x'^, in that rapid spatial variations in electron density, such as those at subsys- 
tem boundaries, can be handled more effectively; LDA-like local correlations are 
also included. Thus the theory has credibility even at small separations outside 



the traditional asymptotic van dcr Waals limit. Wc can therefore hope for a 
consistent description of the forces between subsystems such as molecules or 
metal surfaces, equally reasonable at all separations from the chemical bonding 
regime through to the nonretarded van der Waals limit. 



1. INTRODUCTION 

Jay Mahanty is well known for his pioneering work on van der Waals ( "vdW" 
or "dispersion") forces [1-7], and the present paper is presented as a tribute in 
the form of a possible extension of his methods. Perhaps the most familiar exam- 
ple of a dispersion force is the attractive interaction between a widely separated 
pair of neutral atoms, which is proportional to \/R^ in the nonretarded regime 
and to 1/W in the retarded regime, where R is the interatomic separation. 
More generally, dispersion forces become important in the interaction between 
electrically neutral subsystems at separations larger than a few Angstroms so 
that electron transport between the subsystems is negligible. The present paper 
attempts to generalise existing theories of the van der Waals interaction to cases 
where the separation is not necessarily large compared with atomic dimensions, 
so that that the interaction between the subsystems cannot be treated perturba- 
tively. In this regime one has to consider Angstrom-scale details of the electron 
distribution at the edges of the subsystems. For a description of electronic be- 
haviour on such short spatial scales one cannot rely upon the hydrodynamic 
arguments which have recently been popular in the van der Waals context, so a 
more microscopic theory is required. An additional motivation for the present 
paper is the suggestion [4,8-11] that van der Waals phenomena can be important 
even in large connected regions of high electronic density {e.g. bulk metals) for 
which van der Waals phenomena have not traditionally been considered. 

The van der Waals energy is relevant in the interaction between neutral 
systems in the domain where the separation is greater than a few Angstroms, 
but less than the wavelength of light at the dominant fluctuation frequencies so 
that electromagnetic retardation effects [1,3,5] can be ignored. The vdW inter- 
action can be derived, for example, by second-order perturbation theory in the 
coulomb interaction between the two species [1], or by an RPA-like fluctuation- 
dissipation theorem argument involving the electric polarisability, which will 
be discussed further below. On the other hand, for quantum systems in close 
contac;t. the simplest reasonably successful general theory describing the Pauli 
repulsion and bonding forces is the Local Density Functional theory of Kohn 
and Sham, together with improvements by a number of workers [12,13]. The 
local density functional expression for the exchange-correlation energy misses 
the van der Waals interaction, however, and RPA calculations are not reliable 
for bonding situations nor for extended systems at metallic densities. There 
is therefore a need for a theory which works in all regimes. Ashcroft and co- 
workers [8-11] have recently suggested ways of including van der Waals effects 
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in a density functional framework. In the present paper we suggest a rather 
different quasi-local density formalism which should reproduce the vdW force 
between widely separated species but which remains well-defined as the species 
arc brought together. It makes use of a prior calculation of the dynamic elec- 
tronic susceptibility of the uniform electron gas. 

The new method will now be developed as follows. Firstly, in section 2 the 
RPA groundstate energy formalism is reviewed for the case of an inhomoge- 
neous electronic system: this amounts to a time-dependent Hartree or mean- 
field evaluation of the inhomogeneous dynamic susceptibility x followed by use 
of the fluctuation-dissipation theorem in order to generate approximate pair 
correlations. A coupling-constant integration is also involved. Then in section 
3 it is shown that the vdW interaction between separated neutral systems is 
reproduced by the RPA formulation described in Section 2. In Section 4 it is 
pointed out that the essential element in producing the vdW energy term is not 
the nonlocality of the bare susceptibility (which merely refiects independent 
electron transport, a phenomenon not involved in the vdW interaction). Rather, 
the essence of the vdW interaction stems from nonlocality of the selfconsistent 
Coulomb screening (Hartree) integral equation, which allows distant density 
fluctuations to interact. It is therefore proposed to make a suitable local or 
quasi- local approximation for followed by solution of the nonlocal Coulomb 
screening integral equation to produce an interacting susceptibility x- Finally 
one uses the fliictuation-dissipation theorem and a coupling-constant integra- 
tion. This approach leads to an explicit though fairly cumbersome procedure 
for obtaining a total energy including the vdW contribution, starting from a 
given static density n(r). This expression is merely an approximation to the 
full RPA energy and does not include the local field corrections inherent in, for 
example, the Kohn-Sham exchange-correlation functional. To remedy this omis- 
sion we propose in section 5 a way to include LDA-like local field corrections 
alongside the RPA effects already introduced. A summary is given in Section 6. 



2. SUMMARY OF THE RPA GROUNDSTATE ENERGY FORMALISM 

Since it is being asserted that the vdW interaction is inherent in the inho- 
mogeneous RPA correlation energy formalism, we begin by stating exactly what 
is meant here by this formalism. To calculate the RPA linear response of an 
interacting system one first solves the selfconsistent Hartree groundstate equa- 
tions [14] to yield a selfconsistent potential V°'^{r,s), plus a set of one-particle 
orbitals ipi{r,s) and eigenvalues ej. Then from first-order perturbation theory 
one forms [20] the bare susceptibility 

Xir,r,.) = Y^iU - U) e._e,-na; ' 
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where 3 labels the spin projection eigenvalue of a state. Physically, represents 
the linear response, to an external potential 6V{'r) exp{—iu;t), of a system of 
independent electrons moving in the static potential V^^{r) :i.e. 

5n''"^''P{f,t) =exp{-iu;t) j x°{r,r' ,Lo)5V{r')d^r' (2) 

The general form of (2) is dictated by the fact that it is the most general 
possible linear time-invariant connection between the density perturbation Sn 
and the potential 6V which causes it. The fact that x° is spatially nonlocal 
(i.e. connects different "density" and "potential" points r and r') is solely 
due to the motion ( "transport" ) of individual electrons between points r and 
r' : it has nothing to do with electron-electron interactions since x° refers to 
independent electrons. Only in the very simplest theories is strictly local: 
for example static Thomas-Fermi theory assumes slow spatial variations and 
takes X*' = 6^{f— f')dn{r)/d^{r). The derivation of the exact nonlocal form 
(1) for x" amounts to time-dependent perturbation theory for each independent 
electron orbital i, followed by squaring and summing over occupied states i to 
obtain the density perturbation at r. The energy denominator in (1) is the 
standard denominator appearing in the wavefunction perturbation S-ipi, while 
the matrix elements SVij of perturbation theory are reproduced by the spatial 
integral in (2) acting on two of the wavefunctions in (1). Alternatively (1) 
can be obtained from the real-space version of Feynman diagram theory, being 
the retarded version of the open-bubble-diagram polarisability 7r° [16]. For the 
uniform electron gas, (1) is simply the space Fourier transform of the familiar 
bare Lindhard function. 

While hydrodynamic approximations for x° have the correct global or large- 
distance behaviour, they do not obtain the correct f ~ r ' or high— g behaviour 
embodied in (1). Thus hydrodynamic approximations are less able to describe, 
for example, the detailed effects of rapid electron-density falloff at a metal sur- 
face. In due course we will introduce a quasi-local approximation for x° which 
does retain the correct short-ranged behaviour. 

We turn now to the response of Coulomb-interacting electrons to an exter- 
nally imposed potential perturbation. The Time Dependent Hartree Approxi- 
mation [15], otherwise known as the Random Phase Approximation (RPA) [16], 
amounts to the assumption that the electron density responds, via equation (2), 
not to the external potential JF^^*, but rather to a mean- field potential con- 
sisting of SV^^*^ plus the Coulomb field generated by the instantaneous density 
distribution itself. This plausible "screening" argument entails neglect of the 
quantum or thermal fluctuations of the density about its average value. Equiv- 
alently, it neglects the correlations between the positions of the electrons: once 
an electron is discovered at r, its coulomb repulsion reduces the probability of 
finding other electrons nearby, so that the potential it feels due to the other 
electrons is not precisely that generated by the mean density n{r). 

In general, for a sufficiently weak external potential perturbation 
6V^^*{r)exp{—iu;t), an interacting electron system must experience a linear 
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density response of form 



5n(r, t) = cxp(-iwt) j x{r, r', uj)SV''=\f')d^r '. (3) 

The linear coefficient x in tlie general case will be termed the interacting 
susceptibility and can be regarded as the density response at r due to a potential 
perturbation localised at r '. The above mean- field argument then leads to the 
following Coulomb screening equation for x in the RPA: 

x'^"{X,f,f',co) = x°ir,f',co) 
+ J x°(f,fi,a;)Ae2|fi - f2\-^x^''''{\,f2,f',oj)d^ndW (4) 

Here the interaction strength parameter A is included for later convenience 
and should for the present be set to unity. 

The exchange and correlation energy E^c can now be obtained from the 
interacting susceptibility x in the following fashion. For a general inhomoge- 
ncous electronic system, Exc can be expressed exactly in terms of the density- 
fluctuation correlation function [17,18] as 



E^c = \j dX j d^r J d^r' e^\r-r'\-^ {< 5n{f)5h{f') >x -n{r)5^{f - f')] 

(5) 

where 5n{r) = n(r, t) — n{r, t) is the number-density fluctuation operator. 

In the derivation of (5) it was assumed [17,18] that a A-dcpencient external 
potential V{X, r) is supplied to maintain the density n{r) at its fully-interacting 
value while the electron-electron interaction strength is varied from A = to 
A = 1. Without the A integration, equ. (5) is clearly a contribution to the poten- 
tial energy arising from the correlations (both dynamical and exchange-driven) 
between electron positions. The mutual avoidance of electrons also leads to a ki- 
netic energy contribution because of the uncertainty principle, and qualitatively 
it is this feature which necessitates the integration over the coupling-strength 
A. Formally the A integration is derived from a Feynman-theorem argument 
[17,18]. 

We now use a frequency integral to generate the equal-time correlation func- 
tion C{f,f',t = 0) from its time Fourier transform, and then apply the zero- 
temperature generalised fluctuation-dissipation theorem [19], thus obtaining 



< Sn{r)6n{f') >x= (27r)-i / C{X,r,r",Lu)daj 

J —OO 

pOO 

= —{h/7T)Im / (ia;x(A, r, r uj) 

^0 

poo 

= -{n/n) dsx{X,f,f',is). (6) 

^0 
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Here x(A, r, r lu) is the exact linear response of the interacting homogeneous 
system (with pair potential reduced by a factor A) to an external potential 
6V^^*{r} exp{—iu}t). % is in general defined as in equ. (3) except for the presence 
of a reduced interaction strength: 

dnx{r,t) = exp{~iLjt) jx{\r.r',Lo)5V^''\r')cPr'. (7) 

In (6) the frequency contour has been moved by analyticity arguments to lie 
up the imaginary uj axis where % is purely real. Putting equations (5) and (6) 
together we have the following general exact result: 

E^^ = ^J^ dxj d^rd^r' e^\f-r'\-'^ 

xi^-hn-^J^ x{X,f,f",is)ds-n{r}d^{f-f')^ (8) 

This equation allows us to obtain the exchange-correlation energy of an 
arbitrary inhomogeneous electronic system from a knowledge of the interacting 
density-density response x- Essentially the same formula was introduced by 
Harris and Jones [20] in order to investigate the surface energy of a bounded 
jellium metal, the only difference being that in [20] the density was not held 
constant by a A-dependent external potential, so that an extra electrostatic 
energy term was necessary. 

In the case of independent electrons the exact susceptibility is the bare sus- 
ceptibility x°, and then, since there can be no dynamical correlations, equation 
(8) gives the exact exchange energy Ex- The exact correlation energy is then 
found [20] as the difference between the xc energy and the exchange energy: 



Sc = ^y dxjd^'rj d\' ey-r'\-^ p{x{X,r,f',is)-x\r,?',is))ds 

(9) 

In the present case we use the RPA approximation for x- Thus the RPA 
groundstate energy of a general inhomogeneous system is defined as follows. 
First we obtain the bare or independent-electron susceptibility from equation 
(1), for the actual density n(r). Then we obtain the interacting susceptibility x 
within the RPA by numerically solving the screening integral equation (4), for 
each imaginary frequency to = is and for each coupling strength from A = to 
A = 1. Finally the groundstate xc energy is obtained from (8), or equivalently 
the correlation energy is obtained from (9). 

The energy obtained in this fashion may be termed the RPA energy of the 
inhomogeneous many-electron system at hand. This energy is a generalisation, 
for inhomogeneous systems, of the familiar RPA (ring-diagram) groundstate 
energy of the homogeneous electron gas [16]. It is worth stressing here that 
this energy does include an approximate correlation energy, even though, as dis- 
cussed above, the time-dependent Hartree or RPA-screening equations neglect 
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the correlations between electron positions in the presence of an external dis- 
turbance. The resolution of this apparent contradiction lies in the application 
of the fluctuation-dissipation theorem, which relates correlations of electron po- 
sitions in the absence of an external influence, to the response of electrons to 
the presence of a time-varying external influence. The broken symmetry in the 
latter case causes fluctuating density inhomogeneities to occur, and the mean- 
field interactions between these inhomogeneities are the origin of the nontrivial 
correlation function which emerges when the fluctuation-dissipation theorem is 
applied to the RPA linear response function. This nontrivial correlation func- 
tion is in turn the origin of the van der Waals energy which we show below 
is inherent in the RPA treatment of the groundstate energy of inhomogeneous 
electronic systems. 

Note that it is not being claimed that the RPA-scrcencd Coulomb interaction 
between two charges contains the vdW interaction as such: the vdW interaction 
only emerges as a part of the correlation energy after use of the fluctuation- 
dissipation theorem. 



3. VAN DER WAALS INTERACTIONS BETWEEN SEPARATED SYSTEMS 
FROM RPA CORRELATION ENERGY 

In order to show that the van der Waals interaction is inherent in the in- 
homogeneous RPA correlation energy as described above, we now consider the 
RPA description of the situation shown in Figure 1. 




Figure 1. Typical van der Waals situation. The heavy unbroken line repre- 
sents a fictitious surface dividing three-dimensional space into two regions, TZi 
and TZ2. A localised charge- neutral electronic system <Si lies inside TZi well away 
from the boundary, and ^2 is similarly located within TZ2. The van der Waals 
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part of the interaction between Si and S2 is present even if there is an impene- 
trable barrier between the dotted hnes, making electron transport between TZi 
and 7?-2 impossible, so that x°i2 is zero. In the general case the impenetrable 
barrier is of course absent. 

In Figure 1, two charge- neutral systems iSi and 52 are separated by a large 
distance R such that transport of electrons between Si and ^2 is highly improb- 
able. We partition real space into two regions TZi and TZ2 such that TZi contains 
iSi and TZ2 contains 52, with the boundary between 72-1 and TZ2 well separated 
from both Si and <S2. 

We now split the bare susceptibility x°('')^'i'^) into spatially partitioned 
form: 



X^r-, f ', to) = x?,i(r, r ', to) + Xi,2{r, r',u;) + X2,i{r, r", to) + x^^r, r', to) (10) 
where 

^,0 _ fx" if r e Ki and f'e TCi \ ("i i „^ 

- |o otherwise / \'-'-^) 

■sfi _ ( X°(f,r' ,u) if f e TCi and r 'e TC2 "1 /'11K^ 
Al, 2 — "1^0 otherwise / V'-'-") 

and similarly for x'^2,1 and x''2,2 • 

We seek to solve the screening equations in the case that the two neutral 
systems are widely separated, so that x'^i;2 (at the imaginary frequencies we will 
be using) is exponentially small because of a tunneling barrier. In the special 
case that there is an impenetrable barrier separating TZi and TZ2, so that particle 
transport from TZi to 7^2 is impossible (see dotted lines in Fig. 1), then x°i)2 
and x'^2,1 can even be made exactly zero. This special case emphasises the 
fact that the vdW energy we are about to obtain does not depend on electron 
transport from TZi to 'R2, nor on the presence of electron density in the region 
between Si and S2- 

As a reference consider first the case that TZi contains a neutral species Si, 
but 7^2 is empty. (This case corresponds to figure 1 with <S2 removed). Then the 
screening equation (4) for the interacting susceptibility X^^"'i which we denote 
Xi,i in this situation, becomes 

XiA^^r',uj) =Xi,i{r,f',Lj) 

+ j Xi^{r,ri,i^)e^\ri-r2\-\iAr2,r',uj)Sridh2. (12) 

A similar equation holds for the case that T^i is empty, with (1,1) replaced 
by (2,2) throughout. 

Now consider the full RPA screening equation (4) for the van der Waals prob- 
lem, with both polarisable species <Si and <S2 present as in Figure 1. We write 
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the combined interacting susceptibility as the sum of the individual interacting 
susceptibilities plus a correction, 

X = Xi,i+X2,2 + Ax. (13) 

We spatially partition the Coulomb interaction y^(r, r ') = e^|r — r '|^^ and 
also the interacting susceptibility correction Ax, as in equation (11). Then, us- 
ing * to represent spatial convolution, we can write equation (4) as four separate 
equations, one for each of the (1,1), (1,2), (2,1) and (2,2) cases. The first two 
of these equations are 

Xi.i + Axi.i = Xi,i + (Xi.i * + Xi,2 * Vli) * (Xi,i + Axi.i) 

+(x?,i * + x?,2 * ^2';2) * Ax2,i (14) 

Axi,2 = X?,2 + (X?,i * + X?,2 * ^2'l) * Axi,2 

+ (X?,1 * n',2 + X?,2 * y2,2) * (X2,2 + AX2,2). (15) 

Similar equations hold for Ax2,i and Ax2,2- We assume x'^1,2 and x°2,i are 
negligible, corresponding to the vdW limit in which electron transport from Si 
to ^2 is exponentially supressed. Quantities such as x°i,i *^'^ii2 will be treated 
as perturbations. This is reasonable for well-separated subsystems since x*^!,! 
restricts r to lie deep inside region TZi while the index 1,2 on requires r' 
to lie in region 7^2 so that \r — r'\ is large, making small. The interacting 
susceptibility correction Ax is generated by switching on (^^1,2 +V'^2,i ), and 
so all of its partitioned components Axi,i, Axi,2 etc are of first order in the 
perturbation, or smaller. 

Keeping only zeroth and first-order terms in (15) we find 

{I - X?,i * I^m) * Axi,2 = X?,i * * X2,2 + 2ndorder. (16) 

Using 

(l-x;,i*n!i)-'*X?,i=Xi,i (17) 
from equation (12), we can write (16) as 

Axi,2 = Xi,i * V{^2 * X2,2 + {2nrf order in V{^2, ^2,1}- (18a) 
and similarly 

Ax2,i = X2,2 * V^^i * xi,i + {2nrf order in ¥2^^, V^^}- (18b) 

Now expanding equation (14) to lowest nonvanishing order in ^1,2 and 

^2,11 removing zeroth order terms by using (12) in the form 

Xi,i = X°ia +x'^i,i *Xia: and then using (17) and (18b) we obtain 

Axi.i = Xi,i * Vi^2 * Ax2,i + {higher order} 
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= Xl,l * yi,2 * X2,2 * V2^i * Xl,l 



+{3rd order in V^^, V^^}. (19) 

Equations (18) and (19), along with the (2,2) counterpart of (19), summarise 
what happens to the RPA susceptibiUty when we switch on the Coulomb in- 
teraction between two well-separated neutral subsystems. When used with the 
fluctuation-dissipation theorem and a coupling-constant integration as in (9), 
equations (18) and (19) give rise to a term in the correlation energy which is 
precisely the vdW interaction. To see this, note that the correlation energy 
expression (9) calls for the susceptibility x s-t reduced interaction strength A. 
Thus in (18) and (19) (but not in the explicit \r — r'\~^ factor appearing in (9)) 
we replace V\,\ by XVi^i, by XV'^\,2, and similarly for and ^^2,2- 

Note that x" is independent of A and does not depend on (recall that there 
is a single particle potential V{X,r) holding the density constant). Thus, the 
contribution to the correlation energy (9) which depends on ^1,2 and 1^^2,1 is 

Ti f°° 

AEc = -— ds dXA{X,s) (20) 
27r Jo Jo 

where 

A{X,s) = j d^r' V''{f,r')Ax{X,f,f',is) 

= J d^rlV * Ax(A, is)]p^p = TrlV * Ax(A, is)] (21) 

where the trace is defined with respect to integration over r, and has the usual 
cyclic property. Expanding with respect to the spatial partitioning introduced 
earlier we have 

A{X,s) = Tr[V,^^^ * Axl,l(A,^s)] + Tr[V{^2 * Ax2,i(A,is)] 

+Tr[T/2^i * Axi,2(A, is)] + TriVl^ * Ax2,2(A, is)] 

= Ai,2 + ^24 +^2,2 (22) 

Here, from (19) followed by permutation of the last factor to the front, 

Ai^i = Tr[Vl-^ *xi,i(A,is) * AFj^a * X2,2(A, is) * XV^-^ *xi,i(A,is)] 

= X^Tr[xi,i * * Xi,i * V{^2 * X2,2 * ^2^1] (23) 

It is interesting to note that, in order to obtain all terms in (22) which are 
of second order in ¥"1,2, we needed the second-order result in (19) but only the 
first-order result in (18). 
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Now by differentiating the operator identity e(A) ^ * e(A) = I we can show 
in general that 

|.[e(A)-i] = -e(A)-i*|i*e(A)-i (24) 
and applying this to the case ei,i (A) = 7 — Ax°i,i we find 

d d 
^X\A\is) = ^[er,l *X?,i] 

= -er,l * (-X?,i * * * X?,i 

= Xi.i * *Xi,i (25) 

which reproduces the first three operators in (23). Putting (25) into (23) and 
working similarly with ^2,2 we have 



^,1 + ^2,2 = A^Tri— (xi,i) * V^^^ * X2,2 * ^^2^;! + +Xi,i * ^i';2 * ^(X2,2) * v^^] 

= A2 ^Tr [^1^2 * X2,2 * V^^ * xi,i] (26) 
But from (22) and (18b) 

Ai,2 = A Tr[V[\2 * X2,2 * * xi,i] = ^2,1. (27) 
Thus (22) becomes 



A{\ s) = (A'^ + '2X)Tr[Vl^ * X2,2 * V^^ * xi.i] = ^(A'7^r-[ ]). 

Thus the A integration in (20) can be done analytically, giving 

AE" = -n{2-K)-'^ j d^r d^r' d^rid^r2e^\r - r'r^e^ln - r2\~'^ 

x/ dsxi,i{'r,ri,is)x2,2{r2,f',is). (28) 
Jo 

Here the convolutions have now been written out in full, and the suscepti- 
bilities are evaluated at the full interaction strength A = 1. Equation (28) is the 
expression obtained via perturbation theory by Zaremba and Kohn [21] for the 
van der Waals energy between two arbitrary polarisable systems. This shows 
that the van der Waals interaction is contained naturally within the standard 
inhomogeneous RPA groundstate energy formalism. Of course, although the 
form of (28) is correct, the vdW energy produced by the present RPA formal- 
ism will involve the RPA susceptibilities xia and X2,2 of the two separated 
systems, rather than the exact susceptibilities. 
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As an example of the use of equation (28) consider two point polarisable 
dipoles located at foi and ro2, with scalar polarisabilities ai(u)) and 02(0;). Un- 
der the application of an external potential 

$(r) cxp{—iujt) = —e^^SV^^^lr) exp(— iwf), the dipole moment at tqi is propor- 
tional to the local value of electric field: 

pi = ai(a;)(-V$(foi)) = a(a;)e-i W«^*(foi). (29) 

Since a dipole moment pi is produced by moving an electron through a 
displacement x = —e~^pi, this dipole moment can be represented in the weak- 
field regime by a delta-function-derivative electron number density perturbation: 

Sni{r,t) = 6^{f~ [foi + f/2]) - 6^if- [fbi - x/2]) 

~ -f.V(5^(f- foi) = +e-Vi-V S^{r- foi). 

= e-2a(a;)V^F'=^*(foi).V^3(f - foi) (30) 

We determine the susceptibility of the isolated system Si by demanding that 

the density perturbation (30) is reproduced by the generalised response equation 
(3). Thus we find for a point polarisable dipole at roi 

Xi,i(f,fi,a;) = -e-^ai{uj)VrS^{f-foi).VfJ^{n-foi). (31) 

Putting (31) (plus the equivalent for X2,2(f2,f')) i^^to (28) we find, with 
Einstein summation convention on indices j and k, 

AEc = —h{2Tr)~^ ds ai{is)a2{is)x 
Jo 

djdijd2kd'k\f- f'r Vl - "^2!^^ |r=f2=foi, fi=f'=fo2 
- - 00 

ds ai{is)a2{is), R = \foi - fo2| (32) 
which is exactly the van der Waals interaction (see equ 1.18 of ref [1]). 




4. QUASI-LOCAL APPROXIMATION FOR THE 
INDEPENDENT-ELECTRON SUSCEPTIBILITY x° 

The standard local density approximation [12] for the exchange-correlation 
energy can be obtained by making a quasi-local approximation for the interact- 
ing susceptibility Xj as follows: 

X(A, f, f to) ~ x''''^{\ r, r', a;) = x""'^(A, n = n{f), |f - r'\,co). (33) 
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Here x""*'''(A, n, [Arj, a;) is the susceptibility of a uniform electron gas of 
number density n. Putting the (obviously spatially unsymmetric) Ansatz (33) 
into (8) one readily obtains 

E^c = J d^rn{f)e^c{n{r)) (34) 

where ^ 

e^c{n) = -N-'^{h/2TT) dX j d^r j d^r' Xe'^\r-r'\-'^ 

x[ x''"''H>^,n,\f-f'\,is)ds-nh-\S^{f-f')] (35) 
Jo 

is the exchange-correlation energy per particle of a uniform electron gas of num- 
ber density n. 

An inspection of (34) shows that the LDF approximation cannot yield the 

van der Waals energy of a pair of separated electron systems iSi and ^2 as in 
section 3 above. Consider for example two density concentrations separated in 
real space by a region of zero density as imposed by an impenetrable barrier 
(see dotted lines in Fig. 1). Then the exchange-correlation energy from cqu. 
(34) is the sum of the energies of the two isolated systems, there being no 
contribution depending on the separation because of the local character of the 
density dependence in (34). 

The essential idea proposed here is that we should use an Ansatz like (33), 
not for the interacting susceptibility X) but for the bare susceptibility x*^- (It 
might be objected that screening is well known to reduce the effective range of 
response functions so that one should only make local approximations for inter- 
acting response: this is certainly true for the response F{f, f') of the electron 
density at r to an external point charge density at r'. However Xj although 
commonly termed the "density-density response" , in fact describes the density 
disturbance at r due to a point potential disturbance at r '. In this context it 
is noteworthy that, via the Schrodinger equation, the electron wavefunctions 
respond in a quasi-local manner to the local values of the potential, and not 
directly to the charge density which creates the potential. As a result, both x 
and x'^ arc more localised than F. x is not more localised than x" in general, 
however, and in certain cases x is less local than x'^ ■ this occurs for example for 
real frequencies near a plasmon excitation [22], or for widely separated charge 
concentrations as in the classic van der Waals problem studied above. Specifi- 
cally, for the system Si + S2 studied in section 3 above, the special case with a 
hard wall between TZi and 7^2 has x°i)2 = (a kind of locality) but Xi)2 7^ (a 
kind of nonlocality: see cqu. (18)). 

The essence of the present argument is then to presume a form of locality for 
X°, allowing its values for an inhomogeneous system to be approximated from a 
knowledge of the bare susceptibility of a uniform electron gas. This avoids the 
need to calculate wavefunctions in the inhomogeneous situation. The nonlocality 
of the RPA-screened susceptibility x is maintained by solving an explicitly non- 
local screening integral equation in real space. 
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The simplest way to attempt this approach, which is merely an approxima- 
tion to RPA, is to set 

X°(r',r-",w) ~X°'""'V.(^0,|r'-^'l.'^)- (36) 

where 

X°'"""(n,^,w) = (27r)"^y" x°^(n,fc,w)exp(i£f)ci^fc (37) 

is the usual bare uniform-clcctron-gas Lindhard function (Fcynman diagram 
bubble integral), Fourier- transformed into real space. The simplest form for the 
average density is an unsymmetrical one, 

na.(r,f') = n(r) (38) 

which does however have the advantage of preserving particle number as will 
be discussed elsewhere. One could also use a hydrodynamic approximation for 
x". This approach has been extensively investigated by Mahanty. Summerside 
and others [1,2,3]. To implement the hydrodynamic approach in full one has 
effectively to solve spatial differential equations to obtain x*^- Simple analytic 
results can then be obtained which are accurate for the form of the van der 
Waals interactions at large separation. The present method aims to work even 
for small separations and so must be valid both for slow and for rapid spatial 
variations, whereas hydrodynamics is expected to be valid only in the limit of 
slow variations. As will be discussed elsewhere, the present type of approxi- 
mation correctly obtains some of the short-ranged, high— g response properties 
missed by hydrodynamics, and also avoids the solution of spatial differential 
equations in obtaining x°- 

Note that equations (38), (36), (4) and (8) [or (9)], solved in that order, 
constitute a path from a chosen trial electron density n{r) to an exchange and 
correlation energy E^c which we have argued will contain the van der Waals 
interaction, unlike the usual LDF prescription. This occurs because the nonlocal 
screening integral equation (4) is retained: the precise details of the nonlocal 
behaviour of the true x° which are lost in the ansatz (36) are, as we have 
argued, unimportant in obtaining the vdW energy. The advantage of the present 
approach over a full RPA xc energy calculation is that one obtains the bare 
susceptibility approximately from that of a uniform electron gas, without the 
need to find the inhomogeneous one-electron wavefunctions required as in equ. 
(1) for a full RPA calculation. 

When the above prescription is applied to the uniform electron gas one 
obtains, after a little algebra, the uniform-gas correlation energy as given by 
diagrammatic perturbation theory in the RPA or "ring-diagram" approximation 
(see for example equation 12.23 of [16]). This serves to emphasise that we are 
dealing with an RPA type of approximation, applied in this case however to an 
inhomogeneous situation. 

A potential problem with the above algorithm in the van der Waals context 
is that the Ansatz (36), (38) somewhat overestimates the bare response x°i)2 
of the density in one isolated electron concentration <Si to a disturbance 5V^^* 
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occurring in S2- For example, if r e iSi and r' e S2, (38) and hence (36) 
is nonzero even when there is a hard wall between iSi and ^2 so that the true 
X°(r, f',u}) = x^i )2 is strictly zero. This difficulty is inevitable when the average 
density Uav used in (38) depends only on n{r). This is probably not a serious 
difficulty, however, since the spurious response x^1j2 ^^"'^ (which affects results 
mainly via the first term on the right-hand-side of equation (15)) is of Friedel 
form and hence is of order \fi — r2\~'^- This spurious term competes with 
the right-hand side of equ. (16) which is however of order \fi — r2|~^ as it 
is formally the field due to a net-charge-neutral distribution, giving a dipole 
potential in leading order. Thus under a local approximation for the spurious 
response induced is negligible for large separations. More sophisticated quasi- 
local approximations, ensuring that x°i)2 is not overestimated, will be discussed 
elsewhere. 



5. INCLUSION OF LDA-LIKE CORRELATIONS 

The above theory is only an approximation to RPA, and RPA is known to 
yield poor exchange and correlation energies for most condensed matter systems, 
even though it does contain the van der Waals interaction in the case of widely- 
separated neutral subsystems. Here we propose a simple approximate way to 
remedy this by including an LDA-like local exchange and correlation term which 
largely avoids overcounting the RPA correlations already included. 

We achieve this by adding a local-field term Fxc to the time-dependent 
Hartree screening equation (4), which then becomes 

Xx = X°+X°*{V' + F.cx)*Xx (39) 

F,cx{r,n = 5^{f- r-")A-Vxc(A-^n(r)) (40) 
where one possible choice for f^c is given for A = 1 by 

fxcin) = (f {nebcin)) I d'n?. (41) 

Here e,j.c is the exchange-correlation energy per particle of a uniform electron 
gas of number density n. The choice (41) for fxc makes (39), for A = 1, identical 
with the defining equation of the Time Dependent Local Density Approximation 
[23, 24] which is widely used [23, 24, 25] for the finite-frequency response of 
many-electron systems. To generalise this to the case of a reduced Coulomb 
interaction Ae^|r — r'\~^ as in equation (40) we have used the following scaling 
argument. In a uniform zero- temperature electron gas with Coulomb coupling 
constant Ae^, the relevant quantum length is a^* = A^^/?,^/me^, the effective 
Bohr radius. The energy scale is H* = X^me'^/Ji? , the effective Hartree unit. 
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Thus, in a zcro-tcmpcraturc electron gas with density n and reduced Coulomb 
interaction XV^, the local field term f^c = (P {nexc) / dn"^ is given by 

e.cx{n) = A2e(A-3n), f^cxin) = A" V.c(A-=^n). (42) 

To complete the energy functional, x(A, r, r', w = is) from equ. (39) is then 
integrated with respect to X,u),r and r' as in equ. (9), to yield Ec- 

This theory, when applied to the homogeneous electron gas, gives the fol- 
lowing expression for the corrrelation energy per electron: 

/•I r 47rg2 /-oo 
= -n(27r)-^n-i / dX / / ds[x{X,n,q,is)-x°''{n,q,is)] (43a) 

Jo J Q Jo 

where 

X(A, n, q, is) = x°^(n, q, is){l - [AirX^e^q-^ + A" Vxc(A-3n)]x°^(n, q, is)]-' 

(43b) 

Here y^^ is the bare Lindhard susceptibility of the uniform electron gas [16] . 

The quantity e",, from equs. (43) is not necessarily the same as the correlation 
part of the uniform-gas energy t^c already used in the definition (41) of the 
static fxc- The basic reason for this is that we have ignored the q dependence 
and LD dependence of f^c, corresponding to the Time Dependent Local Den- 
sity Approximation. In principle if we used the fxc{X,q,uj) which made (43b) 
the exact dynamical susceptibility of the uniform gas, then an equation similar 
to (43a) would of course give the exact e^c for a uniform gas. We could even 
contemplate the use of a frequency-dependent but local approximation for the 
inhomogeneous-gas quantity fxc{X,f,f',co) (see for example the approach of 
Gross, Kohn and Iwamoto [26,27]) in order to carry out the present xc energy 
scheme. The exact fxc of a uniform gas is only approximately known, however, 
unlike the energy Cxc which is known essentially exactly from Monte Carlo calcu- 
lations. Probably a simpler route, if we are to avoid the inconsistency between 
e"a;c (equs. 43) and the input values of exc, is to find a q — and w-independent 
function fxcin) (not equal to d'^ {ncxc) / dn^) such that the known {e.g. Monte- 
Carlo-derivcd) exc{n) of the uniform gas is reproduced by equations (43). This 
fxc{n) function for the uniform gas could be found once and for all by solving 
(43), regarded as a nonlinear integral equation for the function fxc- This is 
then the "best" local frequency-independent fxc for use in the present theory, 
in the sense that the exact uniform gas xc energy will be reproduced by the 
theory. Clearly such a frequency-independent fxc would represent some sort 
of average over finite frequencies, unlike equ (41) which can be shown [28] to 
yield the exact static response within Kohn-Sham LDA theory. In initial tests 
and applications it may be simplest to use the definition (41) directly without 
attempting to impose a consistency condition on equ. (43), however. 



6. SUMMARY 
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Wc have suggested a theory which combines elements of the RPA (highly 
nonlocal explicit Coulomb screening) and of local density functional theory 
(quasilocal approximation for the BARE inhomogeneous electron gas suscep- 
tibility, plus a local field correction). As a result, after application of the 
fluctuation-dissipation and Feynman theorems, it produces the van der Waals 
interaction in a natural fashion. The inputs to the theory are a trial inhomo- 
geneous groundstatc density n(r) and a imiform-gas response quantity fxc{n) 
from equation (41) (or, preferably, from a self-consistent solution of equs. (43) 
with e"a;c = Cxc)- Starting from the trial n(r) we form the approximate bare 
inhomogeneous response from (38) and (36). This then determines an 
inhomogeneous x via numerical solution of the spatially inhomogeneous screen- 
ing integral equation (39). The xc energy then follows from spatial, frequency, 
and coupling-strength integrations as in equs. (8) or (9). Equations (38), (36), 
(39), and (9), in that order, thus constitute a path from n{r) to a correlation 
energy Ec : in this sense we have a true density functional. The energy is ob- 
tained from n{r) without calculating Kohn-Sham-like orbitals. Advantages of 
the theory should include the ability to calculate the force between two neutral 
systems, yielding the van der Waals interaction at large separation and LDA-like 
results at small separation, while remaining well-defined and physically reason- 
able at intermediate separations. While our calculational procedure for the total 
exchange-correlation energy is more complicated than the LDA algorithm, it is 
certainly tractable in quasi-onc-dimensional systems such as a pair of three di- 
mensional jellium metals with a vacuum gap separating their parallel surfaces. 
(For this case the formalism can probably be carried out [20, 25] with the full 
bare response from cqu. (1), thus providing a test of the local Ansatz (36) in the 
van der Waals context). For arbitrary three-dimensional situations the linear 
screening integral equation (39) is probably the time-limiting step, and since it 
can be expressed approximately as a matrix inversion problem it is amenable to 
parallel computing techniques. 

One drawback of the present theory is that, for two small well-separated 
neutral systems, although the basic separation dependence of the Van der 
Waals interaction will be reproduced, the coefficient may be too large. This is 
because the quasilocal electron-gas estimate of the individual polarisabilities is 
likely to be an overestimate, being based data for a uniform gas whose energy 
levels are closely spaced, in contrast to the widely-spaced levels of small finite 
systems. This drawback should not apply to large systems such as a juxtaposed 
pair of metal surfaces. In this latter geometry, moreover, there already exist [29] 
efficient algorithms for generating suitable tranforms of some of the necessary 
g-space uniform-gas quantities. 
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